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Summary 

A numerical method for solving the three-dimensional bound- 
ary layer equations for bodies of arbitrary shape is presented. 
In laminar flows, the application domain extends from incom- 
pressible to hypersonic flows with the assumption of chemical 
equilibrium. For turbulent boundary layers, the application 
domain is limited by the validity of the mixing length model 
used. In order to respect the hyperbolic nature of the equa- 
tions reduced to first order partial derivative terms, the mo- 
mentum equations are discretized along the local streamlines 
using of the osculator tangent plane at each node of the body 
fitted coordinate system. With this original approach, it is 
possible to overcome the use of the generalized coordinates 
and therefore it is not necessary to impose an extra hypoth- 
esis about the regularity of the mesh in which the boundary 
conditions are given. By doing so, it is possible to limit, and 
sometimes to suppress, the pre-treatment of the data coming 
from an inviscid calculation. Although the proposed scheme 
is only semi- implicit, the method remains numerically very 
efficient. 

1 INTRODUCTION 

A great number of three-dimensional boundary layer calcu- 
lation methods have been developed in the last two decades. 
Some of them are presented in the synthetic papers of Smith 33 , 
Cousteix 1 ^ and, more recently, Humphreys and Lindhout . 
Although the amount of work done to solve the Prandtl equa- 
tions is substantial, some difficulties remain when the cross- 
flow direction changes in the calculation domain. As it has 
been shown by Wang 35 and Krause 21 this problem comes 
from the nature of the set of the boundary layer equations 
which imposes a CFL type condition to the discretization 
scheme (Cebeci et al®). To fulfil this condition, at least two 
solutions may be proposed: i) to choose a simple numerical 
scheme as, for example, an explicit upwind discretization of 
the crosswise derivatives ; ii) to use an implicit discretization 
of the crosswise derivatives at the unknown station. 

With the first solution, the advancement of the in- 
tegration at a given station always goes in the same cross- 
wise direction and the changes of the crossflow, which appear 
on bodies at incidence, cannot be completely calculated, as 
shown by Cebeci^ ^ 25 , unless a change of the discretiza- 
tion scheme across the boundary layer thickness is allowed 


(Lindhout-Moek 25 ). In the second case, the calculation effort 
is much more important, and therefore reduces the interest in 
using the Prandtl equations (Patel-Baek 31 , Johnston 2 ^). In 
practice, a third strategy exists to conciliate the respect of the 
CFL condition with the efficiency of the numerical scheme. 
Considering only the finite difference methods, Cebeci^ uses 
the standard ” Keller Box” method everywhere it is possi- 
ble and the ” zig-zag” scheme where the crossflow direction 
changes. In this latest scheme, the crosswise advection terms 
are partly written at the calculation station, and partly at 
the upstream station. To overcome some limitations of 
this method, Cebeci^ ^ 11 proposes the “Characteristic Box 
Scheme” which takes into account the existence of character- 
istic directions in the boundary layer equations to limit the 
streamwise integration step in the region where the crossflow 
changes sign in the boundary layer thickness. This leads to 
an extra iteration step at each calculation station. 

The numerical scheme which is presented in this pa- 
per integrates the Prandtl equations along the local stream- 
lines, which are sub-characteristic lines. By doing so, the in- 
tegration proceeds always in the same direction whatever the 
crossflow direction, and the CFL condition is fulfilled, provid- 
ing that the marching step is small enough. As the diffusion 
terms are expressed at the unknown station, the proposed 
method belongs to the semi-explicit type. 

The main originality of the proposed method comes 
from the choice of the space in which the equations are inte- 
grated. Most methods use generalized coordinates in a body 
fitted coordinate system. This needs the calculation of the 
Chris toff el coefficients which introduces an extra hypothesis 
dealing with the regularity of the mesh, while the boundary 
layer assumptions impose only the regularity of the body sur- 
face. To avoid this extra limitation, the discretization of the 
equations at a given station can be done in the tangent plane 
to the surface at this point instead of the actual surface. To 
respect the metric properties of the surface and express the 
co variant derivatives of the velocity, the tangent plane must 
be provided with a particular metric. This is simply done 
by orthogonally projecting the body fitted coordinate system 
and the velocity field on the tangent plane at the considered 
points. 


2 Boundary layer equations 

Body fitted coordinate system 

To set up the boundary layer equations, it is convenient to use 
a body fitted coordinate system (see, for example, Hirschel- 
Kordulla*®). Let x* be the cartesian coordinates of a surface 
point. This point is known by the two parameters A" 1 and A" 2 . 
With t\ the cartesian base vector, the vectors defined by 

^ = a = 1 ’ 2 l=I > 2 > 3 (*) 

are tangent to the body fitted coordinate system. 

The surface base reference frame is obtained by 
adding the unity vector of perpendicular to a? and aj. The 
reference frame (e?,ej, ej) in the vicinity of the surface is 
built as shown in figure 1. Introducing the thin layer assump- 
tion, the metric elements g^ 3 = e^. e~* become independent of 
the A 3 coordinate. 

The boundary layer equations are obtained by ap- 
plying the Prandtl hypothesis to the Navier-Stokes equations 
written in the curvilinear coordinates (AT 1 , i = 1,2,3). 

For an incompressible laminar flow, the boundary 
layer equations read : 

V,I/’ = 0 m = 1,2,3 (2a) 

pU'V t U* = -V a P 

d ( dU*\ 

+ dw{ fl dX*) a = 1 - 2 ( 2b ) 

The covariant derivatives of the velocity are expressed using 
the Christoffel coefficients: 

dU* 

v ' u ° = + r ^ J (3) 

In the equations 2a and 2b, the pressure field is known. It is, 
for example, the wall pressure given by an in viscid calculation. 
The boundary conditions are the no-slip condition at the wall 
and the velocity components (with a = 1,2) at the outer 
edge of the boundary layer. The latter can be obtained from 
the pressure field by integrating the Euler equations at the 
wall. 

Nature of the set of equations 

From the theory of quasi linear differential equations, the 
boundary layer equations 2a and 2b are parabolic because 

of the diffusion terms. It has been shown by Wang 3 " 2 and 
Krause^* that the particular influence of the advection terms 
could be studied from the charact eristic surfaces of the sub- 
set of equations made of the first order derivatives. They have 
shown that the surfaces made of the straight lines perpendic- 
ular to the wall and the stream surfaces are sub-characteristic 
surfaces. This means that the influence domain of a particular 
station is limited by the two surfaces, formed of perpendicular 
lines to the wall, which are tangent to the two most deviated 
streamlines. 


3 Numerical method 

A great number of calculation methods have been developed 
to integrate the boundary layer equations in direct mode, 
i.e. with a prescribed external velocity field. Some reviews 
of these methods can be found in Smith 33 , Cousteix* 4 and 
Humphreys- Lindhout^. Most of these methods are space- 
marching, with an upstream discretization of the advection 
terms. 

Lindhout-Boer 2 ** made a semi-implicit method in 
which the crosswise derivatives along X 2 are explicitly dis- 
cretized in the upstream direction, the other derivatives be- 
ing written implicitly. This allows a change of the crossflow 
direction to be taken into account very simply. The calcu- 
lation step in the streamwise direction is limited by a CFL 
condition. To avoid this constraint, it is necessary to express 
implicitly the X 7 -derivatives. This can be done simply if the 
dependence domains remain in a given side of the mesh lines 
X 1 in the whole calculation domain (fig. 2a). For such flows, 
for example flows over infinite swept wings, the calculation ad- 
vances everywhere in the same direction along the X 1 lines. 
Jelliti*^ and Barberis^ used this technique. For more complex 
boundary layer flows, such methods do not allow accessibility 
to the domains for which the crossflow does not remain in the 
marching direction along the A" 1 lines. 

Lindhout et al . have developed a technique in 

which the choice of the numerical scheme for the crosswise 
derivatives in the X 7 direction and the marching sense along 
these lines depend on the most deviated streamlines through- 
out the boundary layer at the calculation station. This allows 
a certain optimization of the calculation effort by choosing in 
each region the most suitable discretization. 

Other methods have been considered. An effi- 
cient scheme of the “predictor-corrector” type is used by 
Matsuno 2 ^. Wang**® has proposed a“zig-zag” scheme in or- 
der to take into account the dependence domains for the dis- 
cretization of the velocity along the A 2 direction. These terms 
are written partly at the known upstream station and partly 
at the unknown calculation station, on both sides of the corre- 
sponding A 1 line. The stability of this scheme is discussed bv 
Krause‘S . A similar scheme has been used also by Iyer et al. 
and Cebeci^. This author prefers a modified version of the 
“Keller box scheme”, called the “characteristic box scheme” 
which takes into account the dependence domains by using 
the direction of the local streamline in the discretization for- 
mulation. This leads to an extra iteration step at each station 
and a limitation of the marching step in the A 1 direction**. 

Fully implicit techniques in which the A 2 -derivatives 
are written in the unknown plane A 1 = Cste (fig. 2a) can 
be considered. Patel- Baek 3 * and Tassa et a/. 3 ^ use the alter- 
nated direction procedure to solve the equations in a whole 
plane A 1 = Cste . Johnston 2 ^ prefers to sweep only in the 
A 2 direction, which leads to iterative inversion of tridiagonal 
matrices; the unknown quantities being taken at the previous 
iteration. 



Equations along the local streamlines 

In order to respect the physical dependence domains at each 
point of the boundary layer while keeping a single marching 
direction along the ATMines, the momentum and energy equa- 
tions will be discretized along the local streamlines. This also 
allows the use of a unique scheme in the whole calculation 
domain. 

As is usually done in boundary layer calculations, a 
reference length L{X l y X 2 ) is introduced to adapt the grid 
perpendicular to the wall to the boundary layer thickness. 
With the normal coordinate 77 = X J /L( A M ,X 2 ), the bound- 
ary layer equations along the local streamlines read 
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with 


dh 


lp = u3 ~ T,u 'w 1 = 1,2 (5) 

dX 1 being the step size in the main marching direction, d 3(77) 
is calculated using the metric coefficients 


Mv) = ( 3.7 dXM) 
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i,; = 1)2 (6) 


where dX % is a function of 77 obtained from the definition of 
the local streamline parallel to the wall 


dXi_dXi m 

U l U 2 

Vtp is the total variation of the velocity component IP along 
the streamline 



a) body fitted coordinate system 



Figure 2 : Building of the calculation mesh and velocity field 
in the tangent osculator plane. 


vir = <*x;, i, a = 1,2 ( 8 ) 

Osculator tangent plane 

The use of generalized coordinates introduces an extra hy- 
pothesis concerning the regularity of the body fitted coor- 
dinate system which must be regular enough to allow the 
calculation of the Christoffel coefficients. Moreover, as the 
calculation method is of semi-implicit type, the respect of the 
CFL condition leads to the use of a subgrid for the integra- 
tion in the X'-direction. The calculations can be done more 
rapidly if the equations are written in a cartesian coordinate 
system. Due to the local character of the boundary layer prob- 
lem, confined to the vicinity of the body surface, it is not the 
global cartesian frame used to define the surface which will 
be considered, but a local cartesian frame linked to a mesh 
of the body fitted coordinate system in which the boundary 
conditions axe given. 

To build the osculator tangent plane, it will be as- 
sumed that the Christoffel symbols are defined, in order to 
show that the new approach is identical to the classic one, 
but this assumption is not necessary. 

Let O be the node (X, 1 , X]) of the mesh in which the 
boundary conditions are given. The local reference frame at 
this point is T x , i = 1 , 2 , 3 . To integrate the boundary layer 
equations to the next station (X^ +1 ,XJ), it is necessary to 
represent in a cartesian space the neighbouring nodes with 
respect to point O as well as the velocity vectors (fig 3 ). To 
this end, at the point 0 of the surface (S) is associated a 



Figure 3 : Sub-calculation mesh with respect to the repre- 
sentation of the body fitted coordinate system in the tangent 
plane. 



Figure 4 : Discretization of the momentum and energy equa- 
tions. 


point 0 f of an euclidian space (E). The reference frame ( e' , 
i = 1,2) at this point is such that 



This leads to the equality for the metric elements 

9,])o> = 9ij)o (10) 

It can be noted that if the points 0 and O' are identical, the 
euclidian space (E) which has been built is simply the tangent 
plane to the surface at 0. In order to give to (E) the metric 
properties that represent the vicinity of point 0 of (S), we 
impose 

a. = to,, (") 

This allows to represent the body fitted coordinate system in 
the vicinity of point 0 by a curvilinear coordinate system in 
the tangent plane while respecting the distances to the second 
order (fig. 2-b). For this reason, the tangent plane is called 
osculator plane. With the condition (11), the image P’ in (E) 
of a point P in (S) near the point O is given by 


O'P' — (e?) 0 [dx l + 1 ( r jt) 0 dx-'dx* 


( 12 ) 


After the construction of the mesh in the neighbour- 
hood of 0 in the tangent plane, the representation of the 
velocity field is simply done: the velocity vectors are known, 
for example, by their modulus and directions with respect 
to the lines X 1 on the surface. The directions with respect 
to the curvilinear mesh in the tangent plane are assumed to 
be the same (fig. 2-c). Knowing the geometry of the mesh 
and the velocity at the nodes, the calculation of the covari- 
ant derivatives of the velocity is straightforward. With the 


representation which has been adopted, the precision of this 
calculation is of first order . • ^ .... 

The covariant derivative of a vector is an intrinsic 
quantity which does not depend on the reference mesh. This 
quantity exists even if discontinuities of the slope of the co- 
ordinate lines are present. In this case, the Christoffel coef- 
ficients are not defined and the velocity components are dis- 
continuous. Such a configuration can be dealt with if the 
osculator tangent plane is built without using the Christoffel 
coefficients. It can be shown that the construction which has 
been described is equivalent to the orthogonal projection of 
the body fitted coordinate system, and the velocity field, in 

the tangent plane at a given point. This transformation re- 
spects the lengths and the angles to the second order, which 
allows to express the covariant derivatives to the first order. 


Basic equation 


It has been shown that the integration of the boundary layer 
equations could be done in the tangent plane instead of using 
the generalized coordinates. For this reason, the equations 
can be written in cartesian coordinates. For a compressible 
turbulent boundary layer, equations (2a), (4a) and the energy 
equation become 
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and the equation of the streamline parallel to the wall 

dx 1 _ dx 3 
u 1 u 2 
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(15) 


dx 1 is given by the marching step along the x'-lines, roughly 
in the general direction of the flow. du a is the variation of 
the u“-component of the velocity over the distance ds along a 
streamline. The energy equation (13c) is written for the total 
enthalpy K 
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In these relations, only applicable in a cartesian reference 
frame, u e is the modulus of the external velocity and the fric- 
tion velocity . 

Laminar-turbulent transition 

Longitudinal instability mode 

Two criteria are used to predict the onset of transition. Both 
are based on stability calculations for the self-similar Falkner- 
Skan velocity profiles and on the relation proposed by Mack^ 
to link the total amplification coefficient n of the most unsta- 
ble instability waves, at the point of transition, to the turbu- 
lence level of the external flow 


n T = —2.4 In T u - 8.43 


(20) 


and the effective viscosity coefficient is expressed as follows 

Mt ~ M+TMt ( 17a ) 

where fi is the dynamic viscosity coefficient given by the law 
of Sutherland, fit the eddy viscosity coefficient and 7 the in- 
termittency function which is equal to 0 for laminar flow and 
1 in turbulent boundary layer. In the transition region, 7 
depends on the thickening of the boundary layer represented 
by the ratio of the momentum thickness to the momentum 
thickness at the beginning of the transition region, Q/0t^ ■ 

Since the first objective of this study is the valida- 
tion of the numerical technique, including the discretization 
scheme and the use of the osculator tangent plane, a sim- 
ple turbulence model is used. The model is a direct exten- 
sion of the mixing length formulation commonly used in two- 
dimensional flows^, with the damping function proposed by 
Cebeci ® 


In a first criterion proposed by Arnal et al} , the ve- 
locity profile is characterized by the mean value of the Pohl- 
hausen parameter X7, and n is represented as a function of 

(*u - J and ^3 

R»u ~ = -206 exp (25.7 Air) 

[in (16.8 T u )- 2.77 Air] (21a) 
r* 0?, du. 
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To determine the critical value of the momentum thickness 
0 Ucr corresponding to the point x^ , the calculated value 
of 0n is compared to the corresponding value of B Ucrs given 
by the stability diagrams and represented by the correlation 


Quo = exp - 14. 8^ N, = ^- 


( 22 ) 


As soon as R^ u becomes equal to Rs xler ,i the instability waves 
become amplified and is reached. 

The second longitudinal criterion, proposed by 
Arnal^, is a parametric type method. For a given velocity 
profile, characterized by the shape parameter H % , the local 
amplification coefficient a t corresponding to the frequency F, 
is represented as a function of R*, in the form of two half- 
parabols. This allows a simplified representation of the stabil- 
ity diagrams with a minimum number of parameters. Know- 
ing the evolution of the shape parameter H along an external 
streamline, the total amplification coefficient is calculated and 
equation (20) is used to determine the onset of transition. 

Streamwise instability mode 

To predict this mode of transition particular to three-dimen- 
sional flows, two criteria can be used. The first one is an 
extension made by Coustols^ of a criterion originally pro- 
posed by Beasley^. The transition occurs when the Reynolds 
number based on the streamwise displacement thickness 
becomes larger than a critical value which is a function of the 


longitudinal incompressible shape parameter. More precisely, 
this criterion reads 

( 0.106 

R*' t = 95 - 5arctan ( (g,- 2 , 3) ^ 

2.3 < Hi < 2.7 (23a) 

Rs 7 t = 150 H t < 2.3 (23b) 

With this criterion, the influence of the turbulence level of the 
external flow is not taken into account. 

The second criterion, also developed by Coustols and 
Arnal^ ^ requires a more important numerical eSort and can- 
not be detailed here. At each calculation station, the most 
unstable direction e of the velocity profiles in the vicinity of 
the crossflow direction must be determined. The transition 
occurs when the Reynolds number defined with the displace- 
ment thickness in the e direction becomes larger than a given 
value which is a function of the turbulence level of the external 
flow. The number and location of the inflection points of the 
velocity profile in the e direction are also taken into account 
in order to represent the results of stability calculations for 
three-dimensional boundary layers. 

Numerical scheme 

The momentum and energy equations (13b) (13c) are discreti- 
zed in the tangent plane according to the scheme presented 
in figure 4. At the unknown station Q, the diffusion terms 
are written at 3 points and the advection term is taken be- 
tween the points Rk and Rk is the origin at the upstream 
station of the streamline going through the point rj k . At this 
stage, all the quantities are known. Rk is calculated according 
equation (15) assuming a linear variation of the velocity com- 
ponents at the upstream stations. This discretization scheme 
leads, after linearization, to three tridiagonal matrices which 
can be inverted independently to give the two velocity com- 
ponents u l and u 2 and the total enthalpy h , ,. The scheme is 
stable whatever the location of points J?* may be. In prac- 
tice, the marching step along X 1 is limited in order that Rk 
remains between the two adjacent stations K and M of the 
calculation point (fig.4). This constraint is identical to the 
CFL condition of a semi-explicit scheme. 

To complete the integration, the normal velocity com- 
ponent u 3 is calculated using the continuity equation (13a). 

The x 1 -derivatives are taken between the points L and Q and 
the x 3 -derivatives are deduced from the relation 

(&),<“-*> (24 *> 

with x 1 and x 2 the cartesian coordinates in the tan- 
gent plane defined in figure 3. 

At each station X 1 , the boundary layer parameters 
we calculated for all the points in the X 7 direction. This is 
always done in the direction of the increasing values of X 2 , 


whatever the crossflow direction may be, because the calcula- 
tion at a particular station is independent of the neighbouring 
points. The process is repeated in the subgrid calculation in 
the direction up to the station A7 +1 of the body fitted co- 
ordinate system in which the boundary conditions are given. 
At this point, the change of direction a of the coordinate sys- 
tem must be taken into account. Since it is imposed that the 
calculation subgrid coincides with the station A7 +1 , a does 
not have to be necessarily small. This means that slope dis- 
continuities of the reference mesh can be correctly treated. A 
new osculator tangent plane is calculated at each node X 7 of 
the station Xj +l and the calculation process continues. 


4 APPLICATION TO A PRO- 
LATE SPHEROID 

To illustrate some capabilities of the method to predict com- 
plex three-dimensional boundary layers, we will consider the 
prolate spheroid with an aspect ratio equal to 6 at a 10° inci- 
dence. A number of experimental studies have been devoted 
to this case, in particular at the DLR ^ ^9 30 cho- 

sen incidence, the experimental pressure field remains close 
to the analytical inviscid pressure field. Moreover, the stag- 
nation point is sufficiently close to the nose of the body to 
use the simple body fitted coordinate system made of ellipses 
passing through the two poles and circles included in planes 
perpendicular to the symmetry axis of the body. 

In figure 5-c, the light lines show the inviscid stream- 
lines at the wall and the thickest lines represent the friction 
lines for a fully laminar boundary layer. The friction lines 
converge to form the separatrice line^ ^2 Along it, a strong 

thickening of the boundary layer occurs, leading to the aban- 
don of the corresponding calculation line after XjL = 0 8- 
Figure 5-a shows the wall friction lines obtained by taking 
into account the transition phenomenon. With = 1.6 10 6 
and a turbulence level equal to 1.5 1Q“ 3 , the boundary layer 
remains laminar in the windward side up to the separation 
line, and turbulent in the leeward side. In the latter side, the 
accumulation of the friction lines for XjL > 0.7 can be inter- 
preted as a secondary separation. In figure 5-b are plotted the 
friction lines calculated by Meier ct al ^0 from measure- 
ments of the skin friction. At a 10° incidence, the influence of 
the flow separation on the pressure field remains small which 
explains the good agreement concerning the location of the 
separation line in figures 5-b and 5-a. The comparison of fig- 
ures 5-a and 5-c shows the great influence of the transition 
phenomenon. 

The same results are presented in figure 6 at a higher 
Reynolds number of 7.2 10®. The transition to turbulence oc- 
curs sooner, which leads to the displacement of the separation 
line towards the leeward region and suppresses the secondary 
separation. 

In figure 7 are plotted the longitudinal and stream- 
wise displacement thicknesses 5* and 62 as well as the shape 
parameter. They are compared to experimental results ob- 
tained by Meier et al at XjL = 0.64 and 0.71. The develop- 
ment of separation is characterized by a thickening of &i and 



a) with transition criteria, R* = 1,6 10 6 T u = 1.5 10~ 3 
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Figure 8 : Ellipsoide at 10° incidence. 

6 2 , particularly important at X/L — 0.71. The evolution of 
the longitudinal shape parameter H is mainly sensitive to the 
nature of the boundary layer. To perform the calculation with 
the present method, the analytical inviscid flow field has been 
used as well as the experimental pressure field. The influence 
on the results remains small. The most critical point concerns 
the prediction of the transition. The external turbulence level 
is equal to 1.5 10" 3 , as in the experiments. With the present 
calculation methods all the transition criteria have been set 
active and the first one to predict transition is retained. As it 
can be seen in the evolution of H in figure 7, the location of 
the onset of the transition near the windward plane of symme- 
try is not correctly predicted. This is difficult to explain be- 
cause the transition occurs along this line by amplification of 
the longitudinal instability waves which are calculated by the 
second criterion^. Maybe the use of the linear instability the- 
ory along a symmetry line with a divergent flow from this line 
must be questioned. Figure 7, showing the skin friction evo- 
lution along the windward and leeward lines in the symmetry 


plane, indicates that the calculated transition point is located 
at X/L = 0.85 with a turbulence level equal to 1.5 10~ 3 . This 
turbulence level gives the correct location of the transition line 
in the lee side region of the body. Its experimental value is 
estimated between 1 and 2 10 -3 . By taking the largest value 
of turbulence level, the transition occurs at Xj L — 0.73 on 
the windward symmetry line, instead of 0.65 experimentally, 
but it reaches 0.17 on the upper symmetry line. 

Calculation time 

In the present method, the marching step in the X 1 direction 
is limited by the most deviated streamline at a given section. 
This step is also limited with respect to the boundary layer 
thickness 6. For the prolate spheroid, the marching step was 
limited to be in the range 0.6<5m m and 2 £ mox , the minimum 
and maximum values being taken in every section X 1 . With 
this condition, roughly 1000 calculation steps are needed in 
the X 1 direction. With 26 lines in the azimutal direction (for 
a half-body), this corresponds on a CRAY XMP to 10 s for a 
fully laminar boundary layer and 30 3 with all the transition 
criteria. 

5 CONCLUSION 

The three-dimensional boundary layer calculation method 
which has been presented is of semi- implicit type. The ad- 
vection terms are discretized along the local streamlines. The 
dependence domains are thereby satisfied with a simple nu- 
merical scheme. The counterpart is a limitation on the size 
of the marching integration step. Despite this limitation, the 
efficiency of the method remains good due to the reduced 
amount of calculation at each step. This is partly a conse- 
quence of the use of local cartesian coordinates. The dis- 
cretization of the equations in the osculator tangent plane 
allows the existence of slope discontinuities in the body fit- 
ted coordinate system in which the boundary conditions are 
given. It also often reduces or suppresses the pre- treatment 
phase of the data for a calculation. 

Although the application cases which have been pre- 
sented only deal with the prolate spheroid at incidence in 
incompressible flow, the application range of the code is very 
large. It extends from subsonic to hypersonic flows. 

For turbulent boundary layers, the mixing length 
model which is used up to now is restrictive. The introduction 
of transport equation model is being done. It has also been 
tested that the present method can run in the inverse mode 
with only minor modifications. 
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